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In ecology, prey switching refers to a predator's adaptive change of habitat or diet in response to prey abundance. In 
this paper, we study piecewise-smooth models of predator-prey interactions with a linear trade-off in a predator's prey 
preference. We consider optimally foraging predators and derive a model for a 1 predator-2 prey interaction with a tilted 
switching manifold between the two sides of discontinuous vector fields. We show that the I predator-2 prey system 
undergoes a novel adding-sliding-like (center to two-part periodic orbit; "C2PO") bifurcation in which the prey ratio 
transitions from constant to time-dependent. Further away from the bifurcation point, the period of the oscillating prey 
ratio period doubles, suggesting a possible cascade to chaos. We compare our model predictions with data and demonstrate 
that we successfully capture the periodicity in the ratio between the predator's preferred and alternative prey types in 
data on freshwater plankton. Our study suggests that it is useful to investigate prey ratio as a possible indicator of how 
population dynamics can be influenced by ecosystem diversity. 

i Introduction 

.2 ' 

^ Ciliates are a type of protist (eukaryotic single cells with animal-like behavior) that propel using an undulating movement 
' g enerated by small hair-like protuberances (called cilia) that cover the cell body. They are also planktoni^, as they are 
I i transported primarily by currents. Ciliates occur in aquatic environments and feed on small phytoplankton, so they constitute 
kn important link between the bottom and higher levels of marine and freshwater food webs [IHj ■ 

The seasonal dynamics of plankton predators and prey, in which a spring bloom in phytoplankton is followed by extensive 
zooplankton predation and nutrient depletion, is well-documented and arises from a combination of abiotic controls (e.g., 
irradiance and vertical mixing depth) and biotic factors (e.g., competition, grazing, and predation at different levels of a 
food web) [38] ■ Importantly, ciliates and their algal prey populations also vary at shorter-than-seasonal temporal scales. 
During years when the spring bloom lasts for several weeks (equivalent to 15 to 30 ciliate generations), phytoplankton and 
tiliate biomasses exhibit recurring patterns of increases followed by declines [45]. In addition, laboratory experiments and 
environmental observations both suggest that ciliates exhibit negative selection against certain types of prey, because different 
ciliate species benefit differently depending on the match between their feeding mode and the prey species that are abundant 
in the prey community [32:. Therefore, it has been suggested that the driving forces for the sub-seasonal temporal variability 

^ lie in predator-prey interactions between diverse predator and prey plankton communities, particularly during periods of the 
year in which environmental conditions are relatively settled [46] . 

■ Coexistence of species in a shared environment can arise from ecological trade-offs [21]. Phytoplankton have traits that 
characterize their maximum growth or photosynthesis rate, their size, or their susceptibility to predation; such properties 
affect how well an individual performs as an organism. Additionally, performance in ecological function is constrained by 
limited resources. An organism that invests more energy in predation-defense mechanisms cannot invest as much energy 
in growth or other facets. Consequently, because of physiological and environmental constraints, species differentiate to 
occupy niches [21) . In an environment that changes seasonally, trade-offs give rise to different optimal periods of time for 
different species |21j . For example, an increase in temperature can be manifested as a shift in species composition, as species 
better suited to the new governing conditions become more abundant. Consequently, temporal variation of environment 
and trade-offs result not only in species coexistence amidst limited resources but also in temporal variation exhibited within 
a community [21 . Much research has investigated the relationship between diversity and temporal variability in terms of 
changes in total biomass and/or in species composition within a community [30'. More recently, many studies in community 
ecology have emphasized traits as a starting point for investigations of co-occurring species |31[ 129] . 
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A standard Lotka-Volterra model for one predator and two prey that allows diversity in the prey community predicts 
extinction of the prey type that has a smaller capacity to survive |22j . Therefore, such a framework is inappropriate for 
investigations of ciliate-phytoplankton dynamics, because several predators and prey coexist, and there is a known preference 
towards one of the phytoplankton prey [32] . One way to resolve this discrepancy is to examine adaptive predator behavior 
in response to changes in prey densities. There have been many investigations of prey switching [33) . in which predators 
express preference for more abundant prey. For example, prey switching has been demonstrated to promote coexistence of 
competing prey species [42] and to decrease prey competition due to a shared predator [1]. As a result, prey switching has 
been suggested as a candidate mechanism for coexistence in communities with diverse prey [5]. 

In smooth differential-equation models for population dynamics, prey switching can be represented using a Holling type- 
Ill functional response [T7], in which a sigmoid function gives the increased degree of predation with increasing densities 
of a principal prey. In a 1 predator- 1 prey system, a Holling type-Ill functional response corresponds to a situation in 
which predation is low at low prey densities but saturates quickly at a high value when prey is abundant. Such a functional 
response was observed in a system of protist predators and their yeast prey [15j . A low abundance of yeast resulted in 
increased sedimentation of yeast at the bottom of a glass tube and the walls, so prey were safe from predation, which resulted 
in a decrease of predator concentration when the prey concentration was below a threshold. Predator preference towards 
more abundant prey can be examined explicitly by constructing models with multiple prey in which the densities of the 
different prey are system variables. This has been done when switching is considered to be independent of total prey density 
[35] as well as for density-dependent switching [2]. The latter accounts for the expectation that it is difficult at low prey 
densities for a predator to distinguish which prey is more abundant [2] . Using information on which prey type was consumed 
last, van Leeuwen et al. derived a functional response for a predator switching between two prey types. When one prey 
density is constant, this functional response is similar to a Holling type-Ill response. However, it is similar to a Holling 
type-II functional response when incorporated into a model for population dynamics of 1 predator and 2 prey. More recently, 
van Leeuwen et al. generalized this method to the case of a predator switching between multiple prey types [48] . 

An alternative approach to studying prey switching is to suppose that predators behave as optimal foragers instead of 
explicitly incorporating adaptive predator behavior in response to changing prey densities. According to optimal foraging 
theory^ a predator's choice to switch prey depends on prey abundances and which diet composition maximizes its rate of 
energy intake [40] . Optimal foraging theory is based on constant population densities, and it aims to predict when a predator 
will always feed on the preferred prey. An optimal forager also feeds on the alternative prey if doing so does not decrease 
its rate of energy intake [40]. Kfivan [23] showed that an alternative prey can be part of an optimal forager's diet (i.e., that 
the probability that the predator attacks the alternative prey is p g [0,1]). By contrast, the corresponding 1 predator-2 
prey system in which the predator is not an optimal forager predicts that both prey exist and that the predator becomes 
extinct. It has also been shown that models based on optimal foragers adaptively feeding on two prey types that differ 
in their profitability have a larger parameter space in which all species coexist and reduced population density oscillations 
compared to models with non-adaptive predators [23] [SS] [6] . Because similar results have been obtained when considering 
prey growth that is not limited by prey carrying capacities, optimal foraging has been suggested as a possible mechanistic 
explanation for phenomena such as species coexistence [2^. In addition to the choice of whether to include or exclude the 
alternative prey, similar models of optimal behavior in 1 predator-2 prey systems have been built to study adaptive change 
of habitat in a two-patch environment [21] and adaptive change of activity level [2^ . 

When using the principle of optimal foraging, a 1 predator-2 prey interaction can be modelled as a discontinuous dynamical 
system. Piecewise- smooth dynamical systems are a class of discontinuous dynamical systems [5J [5] in which continuous 
temporal evolution of state variables alternates with abrupt events that can be caused by friction, collisions, impacts, 
switching relays, etc. Piecewise-smooth systems arise traditionally in engineering applications, such as an impact oscillator: 
the dynamical behavior of the oscillator evolves continuously in time above an impact surface, but it is occasionally disrupted 
because of an impact between the oscillator and a rigid surface [5]. Piecewise-smooth dynamical systems have also been 
used in some applications in biology. For example, in population dynamics, they have been used for the prey switching 
discussed above as well as to analyze effects on predator-prey population dynamics if a predator population is not allowed to 
be harvested when it is rare [121 128] . Additionally, piecewise-linear dynamical systems have been used to approximate Hill 
functions in models of gene regulatory networks [TH [7] and to describe and quantify synchronization of cattle [H] . 

In the present paper, we apply optimal foraging theory to derive a piecewise-smooth dynamical system that models 
a predator adaptively feeding on its preferred and an alternative prey. We propose that the predator has two different 
feeding modes of the preferred prey. The predator behavior is adaptive in the sense that it can choose, depending on the 
prey densities, to consume the preferred prey either with a large or with a small preference. For simplicity, we assume that 
consuming a small amount of the preferred prey amounts to not consuming it at all. Consequently, when choosing to consume 
the preferred prey to a small extent, the predator feeds only on the alternative prey. Similarly, when the predator chooses to 
feed on the preferred prey with a large preference, the alternative prey is not being predated at all. We take two ecological 
trade-offs into account. First, we introduce a linear trade-off in the parameter that represents the prey preference. Thus, an 
increase in the growth in the predator population as a result of feeding on the preferred prey comes at a cost of the energy 
obtained from the alternative prey. Second, to compensate for the preference, we assume that the preferred prey has the 
advantage of a higher growth rate than that of the alternative prey. We examine the outcome of the population dynamics 
in such a discontinuous 1 predator-2 prey system as we adjust the steepness of the linear degree of preference trade-off. We 
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thereby discover a previously unknown (at least to our knowledge) bifurcation in piecewise-smooth dynamical systems and 
provide a possible link between predator-prey dynamics and ecological trade-offs. 

The remainder of this paper is organized as follows. In Section 13.11 we derive a piecewise-smooth dynamical system for 
an optimally foraging predator and two prey in the presence of prey preference trade-off. In Section 13.21 we investigate 
this system analytically. We determine the regions in which the system's dynamics evolve along the switching boundary and 
derive the flow at the switching boundary using Filippov's method. We also present analytical expressions for the equilibrium 
point at the switching boundary and for the points at which the two vector fields on the different sides of the boundary are 
tangent to the switching boundary between them. In Section 13.31 we investigate the 1 predator-2 prey system numerically. 
We examine the dynamics as we adjust the slope of the preference trade-off and discover a previously unknown bifurcation 
in piecewise-smooth systems. This ("C2P0") bifurcation describes a transition between a center located entirely on the 
switching boundary and a periodic orbit that evolves partly along the boundary and partly outside of it. As the distance to 
the bifurcation point increases, the periodic orbit experiences a period-doubling, which suggests a possible cascade to chaos. 
In Section |3^ we compare simulations of our system with data on planktonic protozoa-algae dynamics. Finally, we discuss 
and comment on our model assumptions, results, and possible generalizations in Section |4] before concluding in Section [5l 

2 Lake Constance Data 

Lake Constance is a freshwater lake, with a surface area of 536 km^, that is situated on the Austrian- Swiss- German border. It 
has been under scientific study for decades, and hourly records of weather conditions — such as temperature, surface irradiance, 
and wind speed — have been collected since 1979. In addition to such abiotic factors, there are also weekly data for biomass 
of several phytoplankton and zooplankton species [31 . 

In its natural state. Lake Constance would be categorized as a lake with low level of productivity if there were no excess 
nutrients from agricultural and other sources of human population in its catchment area. However, the increased sewage 
water treatment since 1979 has reduced the amount of nutrients that enter Lake Constance, which was categorized as a lake 
with an intermediate level of productivity (as measured in terms of abundance of nutrients such as nitrogen and phosphorus) 
when the data were collected [IB]. Nowadays, Lake Constance has been categorized again as a lake with a low level of 
productivity. In our comparison between model simulations and data in Section 13.41 we will consider the biomass data for 
the algal prey of protozoan predators until 1999 that were reported previously in Ref. [43l|44|. We concentrate on the early 
growth season, when the population dynamics are governed predominantly by predator-prey interactions rather than by 
factors such as population dilution due to increased mixing between the upper and lower water masses or strong predation 
of protists by carnivorous zooplankton population (which develops later in the season) [38l |46] . 

The data that we have received refer to abundance (individuals or cells per ml) and biomass (units of carbon per m^) 
of a species obtained at least once in a sample of a few ml to a liter of water in Lake Constance between March 1979 and 
December 1999. When unavailable, the abundance or the biomass have been calculated from each other using species size 
as a conversion factor. Each sample was collected from a water column with an area of 1 m^ and a depth of 20 m. The full 
data set includes over 23,000 observations of 205 different phytoplankton species. We have used a subset of these data for 
the comparison in Section [3.41 In particular, we have selected phytoplankton species that were noted in [46] as representative 
species for preferred and alternative prey of ciliate predators. 



3 Piecewise-Smooth 1 Predator-2 Prey Model 
3.1 Equations of Motion 

We begin the construction of a model for a predator and its preferred and alternative prey by assuming a linear trade- 
off between the preference towards the preferred and alternative preys. In our framework, prey switching occurs because 
the predator can adjust the extent of consumption of the preferred prey. We assume that there is a trade-off in the prey 
preference, which effectively is a trade-off in how much energy the predator gains from eating the preferred prey instead of 
the alternative prey. An increase in specialization towards the preferred prey comes at a cost of the predator population 
growth obtained from feeding on the alternative prey. For simplicity, we assume that the preference trade-off is linear: 

92 = -CLqqi + bq , (1) 

where 51 > is a nondimensional parameter that represents the extent of preference towards the preferred prey, > is 
the slope of the preference trade-off, 6, > is the intercept of the preference trade-off, and (72 > is the extent of preference 
towards the alternative prey. Assuming a linear predator mortality and a linear functional response between the predator 
growth and prey abundance, we define a fitness function that the predator maximizes by the net per capita growth rate 

1 dz . 
R= -— ^ eqipi+ eq2P2 - m, (2) 
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where z is the density of the predator population, pi is the density of the preferred prey, p2 is the density of the alternative 
prey, e > is the proportion of predation that goes into predator growth, and m > is the predator per capita death rate 
per day. 

We obtain the switching condition that describes when the predator chooses consumption with a large preference towards 
the preferred prey (gi = ^11) or consumption with a small preference towards the preferred prey ((71 = gig) to maximize 
fitness by substituting ([1]) into Q and differentiating R with respect to q\: 



dR 
dqi 



(pi - aqP2) e . 



(3) 



Thus, when ^ > 0, the largest feasible qi-^ maximizes predator fitness; when < 0, the smallest feasible gig maximizes 
predator fitness. For simplicity, we assume that qi^ = 0. This implies that the predator switches to the feeding mode of 
consuming only the alternative prey when predator fitness is maximized by having a small preference for consuming the 
preferred prey. We also assume (again for simplicity) that prey growth is exponential. This yields the following piecewise- 
smooth 1 predator-2 prey model: 



Pi 

P2 

z 
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/- = 



(ri - I3iz)pi 

r2P2 
{eqiPipi - ■m)z 

riPi 
(r2 - I32Z)P2 

(e(?2/32P2 - rn)z 



, if h = pi- aqP2 > 



, if h = pi - aqP2 < 



(4) 



where h = pi — aqP2 ~ determines the switching manifold, ri > r2 > are the respective per capita growth rates of 
the preferred and alternative prey, and /3i and /32 are the respective death rates of the preferred and alternative prey due 
to predation. In our numerical simulations in Section 13.31 we take /3i = /?2. Hence, we assume that the predator exhibits 
adaptive feeding behavior by adjusting preference rather than attack rate to the governing prey densities. We will specify (|4]) 
at the switching manifold h = pi — aqP2 = in Section 13.2.31 Production-to-biomass ratio (where the biomass is the mass 
of all living and dead organic matter and the production represents the increase in the biomass produced by phytoplankton 
organisms) calculated from measurements in Lake Constance (see Section [2]) can be used as an index for phytoplankton 
growth. These data suggest that typical values for phytoplankton per capita growth rate vary approximately from 0.2/day 
to 0.6/day during the course of a year. 

Because of the prey preference, the predator exerts more grazing pressure on the preferred prey than on the alternative 
prey. Such an advantage for the alternative prey can be explained by a difference in the use of limited nutrients. For example, 
the alternative prey might invest resources in building defense mechanisms, such as a hard silicate cover, that is difficult for 
the predator to digest. As a result, the alternative prey has fewer resources left for population growth than the preferred prey 
(which does not have as good a defense against the predator). To compensate for the difference in preference, we assume 
that the growth rate of the preferred prey is greater than that of the alternative prey. 



3.2 Analysis of the 1 Predator-2 Prey Model 

In this section, we review some definitions from piecewise-smooth dynamical systems and then examine their manifestation in 
the 1 predator-2 prey system (jl]). We will validate our analytical results (and also obtain additional insights) using numerical 
simulations in Section [3.31 



3.2.1 Basic Definitions from Piecewise-Smooth Dynamical Systems 

The piecewise-smooth 1 predator-2 prey system with Lotka-Volterra interaction terms in ^ contains a jump in the derivative 
of the population densities across the discontinuity boundary /i = 0. This jump makes Q a Filippov system, so its dynamics 
can evolve along the switching boundary. As shown in Fig. [1] there are three types of such dynamics in piecewise-smooth 
systems. When the switching boundary is attracting from both sides of the discontinuity, the system is said to exhibit 
sliding [panel (a)]. Crossing occurs when trajectories starting from one side of the discontinuity pass across the switching 
boundary without following the sliding field on the boundary [panel (b)]. If both of the vector fields point outward from the 
discontinuity, then the region is defined as a region in which escaping occurs [panel (c)]. 

One determines the boundaries of sliding, crossing, and escaping by computing the points at which there is a tangency 
between the vector fields /_ or /+ and the discontinuity boundary h — Q (see Fig. [2]). Determining the tangencies is crucial 
for studying the behavior of a system at a switching boundary, and it constitutes the first step for studying how dynamics in 
a piecewise-smooth dynamical system differ from that in a smooth dynamical system . 

There are three basic tangencies between a piecewise-smooth vector field and a switching boundary. In a fold, a vector 
field has a quadratic tangency with a switching boundary [see Fig. [2] (a)]. In a cusp, the tangency is cubic [see Fig. [2] (b)]. 
Finally, a two- fold occurs when two folds intersect, and there is a quadratic tangency between a switching manifold and 
each side of a vector field [see Fig. [2] (c)]. In Section [3.2.5( we will determine the points at which one of the vector fields 
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(c) escaping 



Figure 1: The three possible types of dynamics in a piecewise-smooth system close to the switching manifold h = 0: (a) 
sliding along the sliding vector field fs, (b) crossing, and (c) escaping (i.e., unstable sliding). 

in (|4]) has a cubic or quadratic tangency with the switching boundary h = pi — aqP2 — 0. Although a two-fold can arise in 
three-dimensional piecewise-smooth dynamical systems, it does not occur in ^ because the crossing and sliding boundaries 
do not intersect on the switching manifold. We will derive the regions of crossing and sliding in Section fS. 2. 2 1 



Figure 2: Three basic tangencies between a piecewise smooth vector field and a switching manifold occur at (a) a fold (i.e., 
a quadratic tangency), (b) a cusp (i.e., a cubic tangency), and (c) a two- fold (i.e., an intersection of two folds). [We have 
shaded the sliding and escaping regions.] 

Studying the aforementioned tangencies in a piecewise-smooth system makes it possible to describe bifurcations that 
can occur when, for example, a limit cycle or an equilibrium point intersects a tangency point on a switching boundary. 
Such bifurcations are examples of discontinuity-induced bifurcations (DIB) [5]. In particular, in a system with three or more 
dimensions, such as the 1 predator-2 prey system ([4]), all generic one-parameter sliding bifurcations occur at either a fold, 
a cusp, or a two- fold. Moreover, because the switching manifold /i = in (U) has codimension 1, all of its one-parameter 
sliding bifurcations can be categorized into 8 different cases (depending on the type of the tangency) [TH| ■ 

If there is a cusp, a trajectory situated entirely in the sliding region has a tangency with the boundary. That is, the sliding 
vector field has a quadratic tangency to the switching boundary. Perturbing the bifurcation parameter from the bifurcation 
point results in the sliding trajectory leaving the switching plane tangentially; because of the cubic tangency with the vector 
field, it returns to the sliding region. This is called an adding- sliding bifurcation because a trajectory that was entirely a 
sliding trajectory becomes a trajectory that consists of two sliding trajectories with a non-sliding segment of trajectory in 
between them [5]. A new sliding segment is thereby added to the trajectory. 

In the next four subsections, we determine the regions of sliding, crossing, and escaping for derive the analytical 
expression for the sliding fiow that is the solution to (j4]) at the boundary h — 0] calculate the equilibrium point and 
associated eigenvalues of the sliding vector field; and compute analytical expressions for points at which there is either a 
quadratic or a cubic tangency between one of the vector fields and the switching boundary. 

3.2.2 Regions of Sliding, Crossing, and Escaping 

The regions in which (jlj undergoes sliding or escaping are located where the component of /+ normal to ft, = has the 
opposite sign to the component of /_ normal to /i = 0. The switching boundary is then either attracting or repelling from 
both sides of the boundary. The condition for sliding or escaping to occur in (jj]) is 



where C denotes the Lie derivative along the fiow / and is defined as = /± • V = x|^^ • Thus, either sliding or 

escaping occurs in ^ for > (ri — r2)^. When < (ri — r2)^, the system in ^ crosses the switching boundary because 
the components of /+ and /_ normal to h = have the same sign. Therefore, trajectories that start from one side of the 
boundary pass through h = without evolving along it. 




Cf^hCf_h = {aqP2f [(ri - ~ z^] < , 



(5) 
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In stable sliding regions, the components normal to /i = of both /+ and /_ point towards the switching manifold and 
trajectories reach sliding motion in finite time. This occurs in (|4]) when 

Cf^ < < £f_ =^ ~z < ri - r2 < z . (6) 

Escaping (i.e., unstable sliding) occurs because the components normal to ft- = of both f+ and /_ point away from h = 0. 
Trajectories that hit these regions are repelled from the switching manifold in finite time. Escaping motion is unattainable 
in simulations in forward time. Escaping occurs in ^ when 

Cf^ > > Cf_ =^ -z > n - r2 > z . (7) 



3.2.3 Sliding Vector Field 

The solution to ([4]) at the discontinuity h = can be expressed using Filippov's differential inclusion [13]. According to 
Filippov's method, the flow of ([4| at pi = aqP2 is determined by a linear convex combination of the two vector fields /_ and 
/+ as follows: 

fs = (1 - + a{x)f+ , (8) 

where 



a{x) 



e [0,1]. 



(9) 



If a = 0, then the sliding flow fs is given by /_, so /i < 0. Similarly, if a = 1, then fs = f+, which implies that h < 0. Thus, 
employing ([5]) and at /i = 0, we see that the dynamics of are governed by the sliding vector field 



A - 2 



{ri +r2- z)pi 
(ri +r2~ z)p2 
eqiPi{ri - r2 + z) + eq2P2{r2 - ri + z) - 2mz 



(10) 



3.2.4 The Equilibrium Point 

The sliding vector field (jlOp has an equilibrium point when fs{pi,p2,z) — 0. The equilibrium of the sliding flow located at 
the switching boundary is called a pseudoequilibrium. For the system Q, there is a nontrivial pseudoequilibrium point at 

_ _ aqm{ri + r2) 



e{qiaqri + 92^2) ' 



V2 = —, ■ r , (11) 

e[qiaqri + q2r2) 



f 1 + f 2 



By evaluating the Jacobian of fs at {pi,P2, ^), we find that its three eigenvalues are Ai = and the complex conjugate pair 
A2.3, which satisfy the characteristic equation 



^2 rn{r2 - ri){qiaq - q2) m(ri +r2) . 

2(gia,ri + (72»-2) 2 

The presence of the eigenvalue Ai = suggests that a periodic orbit in ^ might be born in a novel way. As we will discuss in 
detail in Section [3.3.31 this does indeed turn out to be a previously unknown bifurcation. The eigenvalues A2,3 have negative 
real parts when aq > 92/91, are pure imaginary when Uq = 92/91, and have positive real parts when Uq < 92/91- 

3.2.5 Tangency Points 

At the boundary between sliding and crossing regions, /_ or /+ become tangent to the switching manifold h = pi— aqP2 = 0. 
At a fold, one vector field has a vanishing first Lie derivative Cf^ — (but Cj^ 7^ 0) and a nonvanishing second Lie derivative 
of C^^ = [Cf^)"^ 7^ 0. The other vector field satisfies =0, and the gradient vectors of h and Cf^ must be linearly 
independent of each other [9 . The tangency condition for a fold in (|4]) is then 

Cf_^Cf^^Q =^ z = ±(ri - r2) . (13) 

At a cusp, a vector field has a cubic tangency to a boundary. This occurs when both Cfj_ — and >Cy^ = 0. In addition, 
the sliding vector field has a quadratic tangency to the sliding boundary. Thus, the conditions C'^^^ ^ and 7^ must 
hold. Additionally, the gradient vectors of ft, the first Lie derivative Cf^, and the second Lie derivative are required to 
be linearly independent [5]. For /+ in the 1 predator-2 prey system Q, the quadratic tangency is given by 

C)^ = aqP2 [(ri - z)^ -r\~ z{eqiaqP2 - m)] . (14) 
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Substituting Lf^ = into ([H]) yields a cusp at (pi, 7)2,2) = (m/(e(7i), m/(e(jiag), ri — r2). Similarly, for the condition 
for a quadratic tangency is 

£/_ = aqV2 {{ri - zf - z{m - + r^] . (15) 

Substituting Cf_ = into p3|) yields a cusp at (pi,p2,z) — {aqUi/ {eq2),m/ {eq2),r2 — ri). 

3.3 Numerical Simulations 

We now treat the slope of the preference trade-off as a bifurcation parameter and study the 1 predator-2 prey system @ 
numerically. We are interested in the dynamics of Q as the complex conjugate pair of eigenvalues of the pseudoequilibrium 
crosses the imaginary axis. 

The parameter gives the slope of the tilted switching manifold, and it corresponds biologically to the slope of the 
assumed linear trade-off in the predator's preference for prey. A large corresponds to a situation in which a small increase 
in the predator's preference towards the preferred prey results in a large decrease in that towards the alternative prey. When 

^ 0, a small specialization in consuming the preferred prey requires only a small decrease in how much energy the 
predator gains from eating the alternative prey. At Uq = 0, the preference trade-off no longer decreases (i.e., it is flat), 
which corresponds to a situation in which specialization in one prey has no effect on the preference towards consuming the 
alternative prey. 

We simulate ^ numerically in three different cases: (1) > 92/51, (2) aq = (72/91, and (3) < 92/91- These correspond, 
respectively, to (1) Ai = and the complex conjugate pair A2,3 with negative real parts, (2) Ai = and the complex conjugate 
pair A2,3 with real part 0, and (3) Ai = and the complex conjugate pair A2,3 with positive real parts. To obtain our numerical 
solutions, we use the method developed in Ref. |34| for solving Filippov systems. 

3.3.1 The Equilibrium Point 

From our analytical results, we know that the pseudoequilibrium has a eigenvalue and a complex conjugate pair of eigenvalues 
with negative real part when aq > 92/91- Accordingly, numerical simulations of (U) in this parameter regime have trajectories 
that converge to the pseudoequilibrium in (fTTj) . See Fig. |3]for an example phase portrait. 
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Figure 3: Phase portrait for the 1 predator-2 prey model ^ for parameter values aq = 4, 91 = 1, 92 = 0.5, ri — 1.3, 
r2 = 0.26, e = 0.25, m = 0.14, and f3i — (^2 — 1. The system converges to the pseudoequilibrium (black circle) given by (|lll) . 
The predator's diet is composed of both prey types when the system dynamics evolves along the switching boundary h = 
(shaded) according to the sliding vector field (fTOj) . 



3.3.2 Sliding Centers 

At aq — 92/91, the complex conjugate pair of eigenvalues A2,3 have real part 0. As the sliding vector field is neither attracting 
nor repelling in the linearized dynamics, every entirely sliding orbit is a periodic orbit that surrounds the pseudoequilibrium 
(see Fig. |4|). Thus, the amplitude of a periodic orbit depends on the point at which a trajectory intersects the switching 
surface. In addition to these sliding periodic orbits, there exist periodic orbits that cross the boundary between sliding and 
crossing regions, but which initially are not entirely sliding, that converge slowly (with dimensional simulation times on the 
scale of 10'^ days and greater) to an entirely sliding periodic orbit that is tangent to the sliding-crossing boundary (see Fig.H]). 
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Figure 4: Phase portrait for for parameter values Uq = 0.5, qi = 1, q2 = 0.5, ri = 1.3, r2 = 0.26, e = 0.25, m = 0.14, and 
/?! = /32 = 1- The pseudoequihbrium (black circle) is surrounded by periodic orbits that evolve along the switching manifold 
h ~ m the stable sliding region (shaded). Periodic orbits that reach the boundary between sliding and crossing eventually 
converge to the periodic orbit that is tangent to the sliding boundary (blue trajectory). 



3.3.3 Adding- Sliding Periodic Orbit 

For Gq < q2/qi, the pseudoequihbrium (jlip of the sliding flow is repelling because the complex conjugate pair of eigenvalues 
A2,3 have negative real part. There is also a periodic orbit. From our analytical calculations, we know that the vector field 

has a cubic tangency with the switching boundary at the cusp at {pi,p2,z) = {m/ {eqi),ml {eqiaq),ri — r2). Because of 
the cusp, the local flow near the tangency forces the trajectory of the periodic orbit to first leave the switching boundary 
tangentially and then to return to it. In doing this, the periodic orbit acquires a non-sliding segment before returning to the 
switching manifold ft, = 0. Sec Fig. [5] for an example phase portrait. 

Adding-sliding periodic orbits can be detected using the sliding condition, as there are two separate pieces of sliding 
trajectories when pi/{aqP2) = 1. Between the sliding pieces, there is a segment of non-sliding trajectory when pi/{aqP2) > 1. 
Using numerical simulations of we examine how the amplitude of the adding-sliding periodic orbit measured from the 
pseudoequihbrium {pi,p2, z) in (|11|) scales with the distance to the bifurcation point a^crit = 92/91- For Oqcrit — a? < 0.1, we 
record the amplitude as the difference in the maximum and minimum values oi H — H, G — G, and z — z, where H ^ pi~ a.qP2 , 
H = pi— aqP2, G = ttqPi +P2, and G = UqPi +P2- Based on visual inspection, the scaling near a^crit appears to be linear (see 
Fig. El), and we note for comparison that a linear scaling relation is known to arise for the "generalized Hopf bifurcation" for 
piecewise-smooth dynamical systems discussed in Ref. [37] ■ this context, such a generalized Hopf bifurcation refers to a 
periodic orbit that is born when an equilibrium point of a planar, piecewise-smooth, continuous system crosses the switching 
boundary. Here, by a "continuous system," we mean a piecewise-smooth system in which the trajectories always cross the 
switching boundary without evolving along it. In that generalized Hopf bifurcation, a complex conjugate pair of eigenvalues 
of a piecewise linear system (obtained by linearizing a piecewise smooth system) has a negative real part on one side of 
the switching boundary, real part on the switching boundary, and a positive real part on the other side of the switching 
boundary [37] . 

The above bifurcation in (|4]), in which an adding-sliding periodic orbit arises from a center, has not previously been 
studied (to our knowledge). Because a center transitions to such a two-part periodic orbit, we call this a center to two- 
part periodic orbit ("C2P0") bifurcation. In the C2P0 bifurcation, adding-sliding periodic orbits are born via a two-event 
sequence. First, there is a bifurcation reminiscent of the standard Hopf bifurcation from smooth dynamical systems. This 
arises via the behavior of the eigenvalues of the sliding vector field fs, as the pseudoequihbrium changes from attracting to 
repelling and a periodic orbit appears. Second, as the distance to the bifurcation point increases, the periodic orbit grows 
and an adding-sliding bifurcation ensues [5]. 

3.3.4 Period Doubling 

We show a bifurcation diagram for (|4]) by detecting the local maxima of the quantity pi/{aqP2) > 1 when Oq and bq q2- 
The period- 1 adding-sliding periodic orbit that emerges when aq < 92 /qi period-doubles as we decrease Oq away from the 
bifurcation point. As we illustrate in Fig. [T] this suggests that there is a cascade to chaos as — )• 0. From a biological 
perspective, — > corresponds to the situation in which there is little decrease in the preference towards the alternative 
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Figure 5: Phase portrait for (g]) for a, = 0.4, qi = 1, = 0.5, n = 1.3, r2 = 0.26, e = 0.25, m = 0.14, and I3i = ^2 = 1- 
The system has a periodic orbit that leaves the switching manifold (the stable sliding region is shaded) and returns to it 
because of a cubic tangency between /+ and the switching manifold. Most of the time, the predator's diet is composed of 
both prey types because the system's dynamics evolves along the switching boundary according to the sliding vector field 
(|TOl) . With these parameter values, the adding-sliding periodic orbit has a period of approximately 20 days. In a single cycle, 
the alternative prey is excluded from the diet and Pi/ {aqP2) > 1 for approximately 5.5 days. 



prey if the predator has an increase in specialization towards the preferred prey. 

3.4 Comparison of 1 Predator-2 Prey Model Simulations and Planktonic Protozoa- Algae 
Data 

In this section, we compare the prey ratio of an adding-sliding period- 1 orbit from the 1 predator-2 prey model Q with data 
on cryptophyte and diatom prey collected from Lake Constance in spring. 

Previous studies of the Lake Constance data set suggest that ciliates are the most abundant herbivorous zooplankton 
group in spring, whereas cryptophytes and diatoms are the dominant species groups in the phytoplankton community |45) . 
We use 6 different cryptophyte species as the group of preferred prey and 3 different diatoms species to represent ciliates' 
alternative prey. We include data for phytoplankton whose cell size is sufficiently small in comparison to the size of ciliate 
predators and which dominate the algal community in Lake Constance in spring [45^. For both species groups, we use linear 
interpolation to obtain intermediate biomass values for each day of the year from approximately bi-weekly measurements. 
We then calculate the mean of the 20 interpolated yearly data between 1979 and 1999. In Fig. [51 we compare these data 
with the prey ratio that we obtain from simulations of (j4|). Although our model does not capture the increasing trend, it 
successfully reproduces the periodic pattern early in the growing season when predator-prey interaction can be argued to 
govern the protist-algae dynamics more than physical driving forces in water masses that are rich in nutrients (39j . 

4 Discussion 

We have combined adaptive predator behavior and ecological trade-offs to model the dynamics of a predator feeding on a 
preferred and alternative prey as a piecewise-smooth dynamical system. Our model describes an optimal forager that can 
adaptively change its diet depending on the abundances of its preferred and alternative prey. We assumed a linear trade-off 
in prey preference and analyzed the dynamics of the system as we adjusted the slope of the trade-off. To compensate for the 
preference, we assumed that the preferred prey has a larger growth rate than the alternative prey. 

Our model predicts a steady-state for the predator and prey populations as long as the trade-off in prey preference is 
sufficiently steep. In other words, a steady state arises when a small increase in specialization towards the preferred prey 
would result in a large decrease in how much energy the predator could gain from the alternative prey. As we decrease the 
slope of the preference trade-off, a periodic orbit appears and period-doubles as the slope approaches 0. After the bifurcation, 
the population densities oscillate and the prey ratio is no longer constant. From a biological perspective, a mild trade-off 
suggests that a predator with a small increase in the energy gained from the preferred prey would experience only a small 
decrease in the energy gained from the alternative prey. We considered a linear preference trade-off as a first step in studying 
the effect of such a trade-off in a predator-prey interaction in which prey switching occurs. Although several studies and 
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Figure 6: Amplitude oi H — H (circles), G — G (squares), and z — z (asterisks) of (jlj [for parameter values qi = 1, (72 = 0.5, 
ri = 1.3, r2 = 0.26, e = 0.25, m = 0.14, and /3i = /32 = 1] versus the distance a^crit — a? < 0.1 from the bifurcation point 
ftgcrit = Q2/qi- We illustrate the possible scaling relations for a{H — H) (solid curves), a{G — G) (dashed curves), and a{z — z) 
(dash-dotted curves) for three different values: a = 1 (blue), a — 0.5 (red), and a = 1.5 (black). 



observations support the existence of trade-offs [211 [29], it is not clear what kind of functional form they take. A previous 
model for population dynamics and prey switching studied a convex (i.e., concave down) trade-off between the attack rate 
and the degree of specialization exhibited by a predator [2] . Concave relationships have been formulated for the consumption 
of two prey species and the energy obtained from two prey species in a graphical approach for solving the problem of optimal 
diet [36] ■ To develop a more complete picture of the effects of a preference trade-off on populations dynamics, it would be 
useful to consider generalizations of our model with nonlinear trade-off functions. 

In prey switching, a predator can adopt either a large or a small consumption of the preferred prey type according to 
prey abundances. For simplicity, we assumed that the smallest feasible consumption is 0. Traditionally, prey-switching 
models based on optimal foraging theory assume that a predator either includes or excludes the alternative prey and that 
the preferred prey is always part of the diet [40] . Despite our simplifying assumption, our model simulations always converge 
to an attractor in which the preferred prey is not excluded from the diet. We have also constructed an alternative model 
(in addition to that on which this paper focuses) to describe a predator that eats the preferred prey until it becomes rare, 
at which point it also eats the alternative prey. In the presence of a linear prey preference, this piecewise-smooth dynamical 
system has a pseudoequilibrium with a eigenvalue and a complex conjugate pair of eigenvalues that never crosses the 
imaginary plane as long as the growth rate of the preferred prey is assumed to be larger than that of the alternative prey. 
Our numerical simulations of this model suggest that such a 1 predator-2 prey system in which the preferred prey is always 
part of the diet does not exhibit the bifurcation that we see in the 1 predator-2 prey system ((4|) in which the predator eats 
only one prey type at a time. 

We have also constructed a smooth analog of (dj) using hyperbolic tangent functions. The two systems produce similar 
behavior as long as the slope of the hyperbolic tangent is sufficiently large. In particular, both dynamical systems settle to a 
steady state for a, > 22. and to a periodic orbit for < However, when the slope of the hyperbolic tangent is smaller, the 
smooth system predicts that the predator and the two prey also coexist at steady-state levels (instead of oscillating) when 
aq < We presented the analytical expression for the pseudoequilibrium (1111) of the piecewise-smooth system and we 
calculated the equilibrium point for the smooth system numerically. 

In the present paper, we have interpreted adaptation as flexible feeding behavior of a predator that has two feeding modes 
for consuming its preferred prey, and we have assumed that it can switch between them in response to prey densities. Recent 
studies have shown that rapid adaptation of traits affects ecological interactions and can be observed especially in organisms 
with short lifespans, such as species in plankton communities [4j 120] . Interactions between evolutionary adaptation and 
population dynamics have been studied using both a framework of fast-slow and slow-fast dynamical systems. The former 
approach has given insight into how evolution of traits is seen in the population dynamics by studying the limit in which trait 
evolution occurs on a much faster time scale than predator-prey interaction |10| . In the latter framework, the consequences 
that ecological dynamics can have on the evolution of traits have been studied with a piecewise- smooth dynamical system; 
this assumes that evolution occurs at a much longer time scale than predator-prey interactions [TT] . Formulating a fast-slow 
dynamical-system analog of Q would enable an interesting comparison of piecewise-smooth and fast-slow dynamical systems 
and should be insightful for modelling interactions between evolutionary adaptation and population dynamics. 
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Figure 7: Local maxima in pi/{aqP2) > 1 as — for (jl]) with parameter values ri — 1.3, r2 = 0.26, e = 0.25, m = 0.14, 
and Pi = (^2 = I- 



We chose to use exponential prey growth to simplify analytical calculations. An important generalization is to consider 
logistic prey growth, as nutrient limitation is one of the most important nonliving factors that drive the temporal pattern 
of phytoplankton growth [38]. However, it has been suggested that the importance of nutrient limitation is less pronounced 
than protist grazing at the beginning of a growth season |39] for water masses that are rich in nutrients. Although considering 
logistic prey growth increases the number of parameters in a model (in addition to making analytical calculations significantly 
more challenging), it has two key advantages: (1) it would expand the suitable time window for comparing simulations with 
data in water masses that are rich in nutrients, and (2) its less restrictive assumptions yield a model that is also reasonable 
in principle for water masses with low nutrient levels |39J . 

Our 1 predator-2 prey system ([4]) reproduces the periodicity but does not capture the increasing trend in the scaled prey 
ratio during the early growing season in Lake Constance. We speculate that a generalization that allows both the prey growth 
rates and the preference trade-off slope to be time-dependent might make it possible to also capture the increasing trend. We 
motivate the time-dependent prey growth rate by the fact that the phytoplankton production-to-biomass ratio (results not 
shown) calculated for Lake Constance exhibits a linear increase in the prey growth in spring and a decrease in the autumn 
[14j . To motivate the time-dependent slope of the preference trade-off, we remark that there is some evidence for variation 
of the shape of an ecological trade-off in bacteria in different environmental conditions |19) . 

As a first step towards studying models of adaptive predator behavior and ecological trade-offs, we started from the 
simplest case (i.e., a system with 1 predator and 2 prey). In the model ^ that we have introduced, functional diversity 
is present only in the prey community and it arises as the difference in prey growth rates. Accordingly, we have chosen 
species from a large data set to consider representative prey groups and have left prey competition and predator diversity for 
future work. Ciliates are known to have different modes of predator behavior, and they can be categorized roughly in terms 
of being more or less selective j49]. To illustrate a predator that is more selective, we note that some ciliate species hunt 
as interception feeders that scavenge on food particles and intercept them directly. By contrast, ciliate filter feeders sieve 
suspended food particles and provide an example of less selective predators. Such diversity in the predator community can be 
represented using different preference trade-offs. This could be studied using a a piecewise-smooth dynamical system with a 
dimension higher than 3 (or using a fast-slow dynamical system) and which could have more than one switching manifold. In 
addition, the switching manifolds might intersect with each other. Such a generalization would thus be very interesting (and 
challenging) to study from both biological and mathematical perspectives. We have already motivated the former. From the 
latter perspective, we remark that there does not yet exist a general treatment for bifurcations that arise from intersections 
of switching manifolds when the ambient space has more than 3 dimensions 9 . 

Our model exhibits a novel ("C2P0") bifurcation, in which the dynamics transitions at the bifurcation point between 
convergence to a pseudoequilibrium and periodic adding-sliding oscillations through a center. This emergence of adding-sliding 
orbits differs from the usual way that they emerge in piecewise-smooth dynamical systems [5]. The standard mechanism 
entails the birth of an adding-sliding periodic orbit following a bifurcation in which the eigenvalues of the pseudoequilibrium 
cross the imaginary plane, so that the pseudoequilibrium changes from attracting to repelling and an entirely sliding periodic 
orbit is born. As the amplitude of the periodic orbit (which lies entirely on the switching boundary) grows, the sliding 
periodic orbit eventually becomes tangent to the boundary between sliding and crossing. Finally, this periodic orbit becomes 
a trajectory that has two sliding segments separated by a non-sliding segment [5] . We observed numerically that the amplitude 
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Figure 8: Scaled prey ratio pi/{aqP2) for simulations of (|4]) (dashed curves) for — 0.4 (blue), — 0.2 (red), ri — 1.3, 
r2 = 0.26, e = 0.25, m = 0.14, and Pi — 1^2 = ^ and the mean data (solid curves) over 1979-1999 for the prey ratio pi/{aqP2) 
collected in spring in Lake Constance. The preferred prey pi is composed of data for Cryptomonas ovata, Cryptomonas 
marssonii, Cryptomonas reflexa, Cryptomonas erosa, Rhodomonas lens, and Rhodomonas minuta. The alternative prey 
group p2 is composed of data for small and medium-size Chlamydomonas spp. and Stephanodiscus parvus. These data, 
previously reported in Ref. [43l |44] were kindly given to us by Ursula Gaedle. 



scaling of the adding-sliding periodic orbit emerging from the C2P0 bifurcation appears to be linear in the distance from the 
bifurcation, which is also the case for the generalized Hopf bifurcation in piecewise-smooth dynamical systems in Ref. |37j . 

5 Conclusions 

We combined two ecological concepts — prey switching and trade-offs — and used the framework of piecewise-smooth dynamical 
systems to develop a model of one predator feeding on a preferred and alternative prey. We derived analytical expressions 
for the pseudoequilibrium, its eigenvalues, and the points for tangencies between the two vector fields and the switching 
boundary. We confirmed our analytical results with numerical simulations, and we discovered a novel bifurcation in which 
an adding-sliding periodic orbit is born from a center. Based on numerical simulations close to the bifurcation point, the 
amplitude of the adding-sliding periodic orbit seems to scale linearly with the distance from the bifurcation point. 

Our model introduces a way to link trade-offs with adaptive predator behavior. We compared the results of our simulations 
with data on freshwater plankton, but we remark that similar prey-switching models can also be formulated for any other 
1 predator-2 prey interaction in which it is viable to use models based on low-dimensional differential equations (i.e., large 
population size, well-mixed environment, and the use of community-integrated parameters). We also discussed several 
biologically motivated generalizations of our model. We believe that our current model and these generalizations provide a 
promising direction for examining possible mechanisms for ecological trade-offs in population dynamics. 

Although we have motivated our investigation primarily from a biological perspective, it is also important to stress the 
utility of our model for development of new theoretical understanding in piecewise-smooth dynamical systems. The biological 
and mathematical motivations complement each other very well, and investigating the condition for sliding corresponds to 
studying a scaled prey ratio, and this in turn offers a possible link between ecological trade-offs and population dynamics. 
We believe that our model offers an encouraging example of how combining theoretical and practical perspectives can give 
new insights both for the development of theory of piecewise-smooth dynamical systems as well as for the development of 
models of population dynamics with predictive power. 
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